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A  molecular-level  computational  investigation  is  carried  out  to  determine  the  dynamic  response  and  material  topology  changes 
of  fused  silica  subjected  to  ballistic  impact  by  a  hard  projectile.  The  analysis  was  focused  on  the  investigation  of  specific  aspects 
of  the  dynamic  response  and  of  the  topological  changes  such  as  the  deformation  of  highly  sheared  and  densified  regions,  and 
the  conversion  of  amorphous  fused  silica  to  Si02  crystalline  polymorphs  (in  particular,  a-quartz  and  stishovite).  The  topological 
changes  in  question  were  determined  by  carrying  out  a  postprocessing  atom- coordination  procedure.  This  procedure  suggested 
the  formation  of  stishovite  (and  perhaps  a-quartz)  within  fused  silica  during  ballistic  impact.  To  rationalize  the  findings  obtained, 
the  all-atom  molecular-level  computational  analysis  is  complemented  by  a  series  of  quantum-mechanics  density  functional  theory 
(DFT)  computations.  The  latter  computations  enable  determination  of  the  relative  potential  energies  of  the  fused  silica,  a-quartz 
and  stishovite,  under  ambient  pressure  (i.e.,  under  their  natural  densities)  as  well  as  under  imposed  (as  high  as  50  GPa)  pressures 
(i.e.,  under  higher  densities)  and  shear  strains.  In  addition,  the  transition  states  associated  with  various  fused-silica  devitrification 
processes  were  identified.  The  results  obtained  are  found  to  be  in  good  agreement  with  their  respective  experimental  counterparts. 


1.  Introduction 

The  present  work  deals  with  the  problem  of  formation 
of  crystalline  phases  (specifically  a-quartz  and  stishovite) 
within  a  fused-silica  (ceramic  glass  containing  a  high  purity 
Si02)  target  during  a  ballistic  impact.  Hence,  the  main  aspects 
of  the  present  work  include  (a)  short-  and  intermediate-range 
order  topology  of  fused  silica,  (b)  (crystalline)  polymorphs 
of  Si02,  and  (c)  devitrification  (i.e.,  crystallization  of  glass 
taking  place  under  the  ballistic-impact  loading  conditions). 
A  brief  overview  of  these  aspects  of  the  problem  at  hand  is 
presented  in  the  remainder  of  this  section. 

Short-  and  Intermediate-Range  Order  Topology  of  Fused  Silica. 
Ceramic  glasses  (such  as  fused  silica),  like  metallic  glasses  and 
glassy  polymers,  are  amorphous  materials.  The  molecular- 
level  topology  of  these  materials  involves  entities  such  as  a 
random  network  of  covalently  bonded  atoms,  atomic  free 
volumes,  (network  former)  bridging  and  nonbridging  oxy¬ 
gen  atoms,  and  cations  of  glass-modifier  species.  However, 
despite  the  absence  of  a  crystalline  structure,  the  topology  of 


ceramic  glasses  is  not  completely  random.  Rather,  it  involves 
different  extents  of  short-  and  intermediate -range  order 
spanning  over  a  range  of  length-scales  (from  the  quantum- 
mechanical  to  the  continuum-level).  To  describe  the  struc¬ 
ture  of  ceramic  glasses  as  determined  using  various  exper¬ 
imental  techniques  (e.g.,  neutron-diffraction,  nuclear  mag¬ 
netic  resonance,  and  small  angle  X-ray  scattering  (SAXS)), 
the  so-called  “random  network  model”  [1]  is  typically 
employed.  Such  a  model  represents  an  amorphous  material 
as  a  three-dimensional  linked  network  of  polyhedra.  The 
character  (number  of  facets)  of  the  polyhedra  is  controlled 
by  the  species- specific  coordination  of  the  central  (glass- 
forming)  atom  (cation).  In  the  case  of  silicate-based  glasses 
like  fused  silica  and  soda-lime  glass,  the  polyhedron-center 
atoms  are  all  silicon  and  each  silicon  atom  is  surrounded  by 
four  oxygen  atoms  (while  each  oxygen  atom  is  connected  to 
or  bridges  two  silicon  atoms)  forming  a  Si044-  tetrahedron. 
In  glasses  of  different  formulations,  other  types  of  polyhedra 
may  exist.  Since  silicon  has  a  tendency  to  form  a  continuous 
network  with  (bridging)  oxygen  atoms,  Si02  is  commonly 
referred  to  as  a  “network  former.” 
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(a)  (b) 

Figure  1:  The  atomic  arrangements  within  the  nonprimitive  unit  cells  of  (a)  a-quartz  and  (b)  stishovite. 


Polymorphs  of  Si02.  Previous  investigations  (e.g.,  [2-14]) 
established  that,  under  high-rate  (shockwave  or  ballistic) 
loading  conditions,  fused-silica  targets  can,  at  least  in  the 
vicinity  of  the  impact  region,  experience  transformation  of 
the  amorphous  structure  into  a  crystalline  one.  Examination 
of  the  standard  temperature-pressure  phase  diagram  for  Si02 
reveals  that,  at  room  temperature  and  slightly  above  it, 
Si02  can  appear  in  one  of  the  following  three  polymorphs: 
a-quartz,  coesite,  and  stishovite.  Consequently,  during  a 
ballistic  impact  into  a  fused-silica  target,  these  three  Si02 
polymorphs  are  most  likely  to  form  (provided  such  impact 
produces  a  local  crystallization  of  the  fused  silica).  Since  the 
main  crystallographic  features  of  these  three  polymorphs  can 
be  readily  found  in  the  public-domain  literature,  they  will  not 
be  overviewed  here  but  will  be  used  in  the  computational 
analyses  presented  later  in  this  paper.  The  atomic  arrange¬ 
ments  within  the  nonprimitive  unit  cells  of  (a)  a-quartz  and 
(b)  stishovite  are  given  in  Figures  1(a)  and  1(b),  respectively. 
The  same  is  not  given  for  coesite  since  this  Si02  polymorph 
does  not  generally  form  under  ballistic-impact  conditions. 

Dynamic-Loading-Induced  Crystallization  of  Glass.  A  detailed 
overview  of  the  public-domain  literature  carried  out  as  part 
of  the  present  work  revealed  a  number  of  experimental  and 
computational  investigations  dealing  with  the  mechanical 
response  of  fused  silica  to  dynamic  loading.  Some  of  these 
studies  revealed  the  formation  of  shear  bands  within  other¬ 
wise  amorphous  glass  and  others  established  the  formation  of 
crystalline  phases  (a-quartz  and  stishovite,  but  not  coesite), 
while  still  others  demonstrated  increased  hardness  of  the 
material  surrounding  the  impact  region  without  establishing 
the  topological  cause  for  this  property  change. 

Chakraborty  et  al.  [15]  carried  out  a  series  of  nanoinden¬ 
tation  tests  and  showed  that  as  the  loading  rate  increases,  the 
extent  of  shear  band  formation  in  the  region  surrounding 
the  indentation  decreases,  while  the  hardness  value  increases. 
No  crystal- structure  analysis  was  carried  out  to  determine 


potential  formation  of  any  of  the  crystalline  phases  as  a  result 
of  loading  or  to  provide  a  rationale  for  the  observed  effect  of 
the  loading  rate. 

Tschauner  et  al.  [16]  investigated  formation  of  stishovite 
in  soda-lime  glass  during  57  GPa  shock  loading  experi¬ 
ments  and  the  reversion  of  this  phase  during  subsequent 
release/unloading.  They  demonstrated  that,  upon  loading, 
high-density  non-fully-crystallized  Si02  phase  was  present  in 
the  “shocked”  fused  silica.  Upon  static  loading  to  only  13  GPa, 
this  phase  was  fully  converted  into  the  crystalline  stishovite, 
suggesting  that  the  shock  loading  was  able  to  devitrify  fused 
silica  and  form  crystalline  stishovite. 

Salleo  et  al.  [17]  demonstrated  the  formation  of  a  defective 
form  of  stishovite  in  the  surface  region  of  fused-silica  wafers 
through  irradiation  with  a  high-powered  laser.  The  formation 
of  such  a  phase  and  its  continuous  growth  has  been  found 
to  be  the  main  cause  of  failure  in  optics  used  for  high-power 
photonics. 

Mantisi  et  al.  [18]  carried  out  a  comprehensive  atomic- 
scale  simulation  of  fused  silica  under  combined  pres¬ 
sure/shear  loading  conditions.  The  results  obtained  estab¬ 
lished  permanent/irreversible  densification  of  the  fused- 
silica  test  sample  and  the  change  of  the  silicon  and  oxygen 
coordination  relative  to  that  present  in  the  as-received  fused 
silica.  However,  these  topological  changes  did  not  appear  to 
be  the  result  of  glass  devitrification/crystallization  processes; 
that  is,  the  glass,  while  compacted,  remained  amorphous. 
One  of  the  potential  reasons  for  the  differences  in  the  findings 
reported  in  the  present  paper  and  in  the  work  of  Mantisi 
et  al.  regarding  the  response  of  glass  to  high  pressure  and 
shear  is  the  use  of  different  interatomic  force-field  potentials. 
It  is  well-established  that  relatively  minor  modifications  in 
the  force-field  potentials  may  have  profound  effects  on  the 
outcome  of  the  all-atom  molecular-level  simulations. 

Kubota  et  al.  [19]  used  molecular  dynamics  simulations 
to  infer  the  atomic-scale  structural  changes  in  fused  silica 
induced  by  shock-compression  loading.  The  results  obtained 
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revealed  that  shock-compressive  loading  involving  stress  lev¬ 
els  exceeding  the  Hugoniot  Elastic  Limit  gives  rise  to  dramatic 
changes  in  the  structure  and  topology  of  the  fused- silica 
network  and  densifications  in  excess  of  20%.  Coordination 
analysis  of  the  as- shocked  fused  silica  revealed  the  formation 
of  under-  and  overcoordinated  Si  atoms.  While  undercoordi¬ 
nated  Si-atom  regions  could  be  interpreted  as  shock-induced 
fused-silica  flaws,  overcoordinated  Si-atom  regions  showed 
some  resemblance  to  the  crystalline  stishovite.  In  addition 
to  the  coordination  changes  just  described,  changes  in  glass 
topology  (such  as  increases  in  the  number  of  the  threefold, 
fourfold,  sevenfold,  and  larger  rings)  were  observed. 

Main  Objectives.  The  main  objective  of  the  present  work  is  to 
carry  out  a  series  of  all- atom  molecular-level  computational 
investigations  of  the  ballistic  impact  by  a  hard  projectile  onto 
a  fused-silica  target-plate  and  to  characterize  the  target-plate 
material  response  to  such  impact.  In  addition,  in  order  to 
help  rationalize  some  of  the  findings  regarding  the  topology 
of  the  fused  silica  following  the  impact,  a  series  of  quantum- 
mechanical  density  functional  theory  (DFT)  analyses  will 
be  carried  out.  These  analyses  will  help  reveal  the  relative 
potential  energies  of  the  Si02  amorphous  state  and  two  Si02 
crystalline  polymorphs,  that  is,  a-quartz  and  stishovite,  and 
the  changes  in  these  energies  as  a  function  of  the  extent  of 
material  compression  and  shear. 

Paper  Organization.  Details  regarding  the  all-atom  mole¬ 
cular-level  computational  procedure  used  to  simulate  the 
ballistic  impact  by  a  hard  projectile  onto  a  fused-silica  target- 
plate  are  presented  in  Section  2.  The  quantum-mechanical 
DFT  procedure  used  to  determine  the  relative  potential  ener¬ 
gies  of  the  Si02  amorphous  state  and  the  two  Si02  crystalline 
polymorphs,  as  well  as  the  associated  transition  states  (to  be 
defined  later),  is  provided  in  Section  3.  Key  results  obtained 
in  the  present  work  are  presented  and  discussed  in  Section  4, 
while  a  summary  of  the  main  findings  and  conclusions  is 
provided  in  Section  5. 

2.  Molecular-Level  Analysis  of  Ballistic  Impact 

As  mentioned  earlier,  one  of  the  objectives  of  the  present 
work  is  to  carry  out  a  series  of  all-atom  molecular-level 
computational  analyses  of  the  problem  of  a  ballistic  impact 
by  a  hard  projectile  onto  a  fused-silica  target-plate.  Within 
the  all-atom  molecular-level  computational  methods  and 
tools,  every  atom  and  every  bond  is  explicitly  accounted 
for  and  molecular  mechanics  and  dynamics  algorithms  are 
used  to  quantify  the  state  and  behavior  of  the  material  under 
investigation.  All-atom  molecular-level  simulation  problems 
typically  require  the  specification  of  the  following:  (a)  a 
molecular-level  computational  model  consisting  of  atoms, 
ions,  functional  groups,  and/or  molecules;  (b)  a  set  of 
force-field  functions  (mathematical  expressions  describing 
bonding  and  nonbonding  interactions  between  the  model 
constituents,  e.g.,  atoms  and  ions);  and  (c)  computational 
method(s)  to  be  used  in  the  simulation.  More  details  of  these 
three  aspects  of  the  molecular-level  modeling  and  simulation 
of  fused  silica  are  provided  below. 


2.1.  Computational  Model.  The  computational  model  used 
in  this  portion  of  the  work  consists  of  two  distinct  subdo¬ 
mains:  (a)  the  projectile  subdomain  and  (b)  the  target-plate 
subdomain.  The  two  subdomains  are  shown  and  labeled  in 
Figure  2(a). 

The  projectile  subdomain  is  of  a  right-circular  solid  cylin¬ 
drical  geometry  (height  over  diameter  ratio  =  1.0,  axis  aligned 
with  the  ^-direction).  While,  at  least,  the  core  of  ballistic 
projectiles  is  typically  made  of  hard  and  heavy  metallic 
materials  (e.g.,  tungsten),  metallic  materials  could  not  be  used 
in  the  construction  of  the  projectile  in  the  present  work.  The 
reason  for  this  is  the  absence  of  metallic  force-field  functions 
(in  the  pure  metallic  environment)  within  the  force-field 
function  database  used  in  the  present  work.  Consequently, 
the  projectile  was  made  of  another  hard  material,  diamond. 
Typically,  the  projectile  subdomain  contained  150  covalently 
bonded  carbon  atoms  forming  a  perfect  single-crystalline 
diamond  structure. 

As  far  as  the  target-plate  subdomain  is  concerned,  it  is  of 
a  rectangular  parallelepiped  (plate-like)  shape  and  is  made  of 
fused  silica.  At  the  molecular  level,  fused  silica  is  modeled 
as  a  discrete-particle  based  material  consisting  of  silicon  (Si) 
and  oxygen  (O)  atoms  mutually  bonded  via  a  single  covalent 
bond  and  forming  a  connected,  nonstructured/amorphous 
network  of  silica  (Si044-)  tetrahedra.  While  fused  silica  is 
an  amorphous  material  and  does  not  possess  any  long-range 
regularity  in  its  atomic/molecular  structure,  modeling  of  bulk 
behavior  of  fused  silica  is  typically  done  at  the  molecular 
level  by  assuming  the  existence  of  a  larger  (amorphous) 
unit  cell.  Repetition  of  this  cell  in  the  three  orthogonal 
directions  (the  process  also  known  as  application  of  the 
“periodic  boundary  conditions”)  results  in  the  formation 
of  an  infinitely  large  bulk-type  material.  This  procedure 
was  adopted  in  the  present  work.  The  parallelepiped-shaped 
target-plate  computational  subdomain  used  in  the  present 
(ballistic-impact)  analysis  contained  9600  particles  (3200 
Si  atoms  and  6400  O  atoms).  The  edge-lengths  of  the 
computational  cell  were  initially  set  to  a  -  b  -  7.9  nm,  c  ~ 
2.9  nm  (approximately),  yielding  a  fused-silica  initial  nomi¬ 
nal  density  of  2.19  g/cm3.  The  three  edges  (a,  b,  and  c)  of  the 
cell  were  aligned,  respectively,  with  the  three  coordinate  axes 
(x,  y,  and  z)  (with  the  target-plate  thickness  aligned  with  the 
-^-direction). 

To  create  the  ambient  temperature/pressure  “equilibrium” 
atomic  configuration  within  the  computational  cell,  the 
following  procedure  was  implemented  within  the  Visualizer 
[20]  program  from  Accelrys. 

(a)  A  starting  computational  cell  was  first  constructed 
by  stacking  the  appropriate  number  of  a-quartz  unit  cells 
in  the  three  orthogonal  directions.  This  was  followed  by  the 
affine  distortion  of  the  computational  cell  in  order  to  obtain 
the  correct  mass  density  of  the  fused-silica  amorphous  state 
(obtained  using  the  following  procedure). 

(b)  To  convert  the  crystalline  material  into  an  amor¬ 
phous  one,  a  stochastic  bond-switching  algorithm  was  then 
implemented  [21]  using  a  Monte  Carlo  computational  pro¬ 
cedure.  Within  this  algorithm,  two  neighboring  Si-O  pairs 
were  randomly  selected  from  the  computational  cell  and 
the  potential  energy  change  A E  resulting  from  the  Si-O 
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Figure  2:  (a)  The  computational  model  used  in  this  portion  of  the  work,  consisting  of  two  distinct  subdomains  for  the  projectile  and  the 
target-plate,  and  (b)  close-up  of  the  resulting  fused-silica  molecular-level  random  network  topology.  Larger  ball  sizes  are  used  in  this  figure 
to  highlight  a  pair  of  Si044~  tetrahedra  sharing  a  common  oxygen  atom. 


bond  switching  computed.  In  the  AE  <  0  case,  the  bond 
switching  in  question  was  accepted  without  any  additional 
conditions.  On  the  contrary,  in  the  AE  >  0  case,  a  Boltzmann 
probability  factor  PB  =  exp[-AE/(3NkT/2)]  (where  N  is 
the  number  of  atoms  within  the  computational  cell,  k  is  the 
Boltzmann  constant,  and  T  is  the  absolute  temperature)  was 
first  calculated  and  compared  with  a  random  number  RN 
drawn  from  a  (0, 1)  uniform  distribution  function.  The  bond 
switching  in  question  was  then  adopted  only  if  PB  >  RN. 

(c)  The  resulting  structure  was  then  subjected  to  a 
carefully  devised  set  of  NVT  (where  N  (=9600)  is  the  (fixed) 
number  of  atoms  within  the  computational  cell,  V  is  the 
computational  cell  volume  (also  fixed),  and  T  is  a  fixed 
temperature)  molecular  dynamics  simulations.  Specifically, 
the  NVT  simulations  were  started  at  a  temperature  of  5300  K 
and  carried  out  in  such  a  way  that  the  temperature  was 
controlled  using  the  following  “  simulated- annealing”  scheme: 
(i)  a  particle-velocity  scaling  algorithm  was  applied  every 
time  step  for  the  first  6000  steps.  This  enforced  strict  control 
of  the  temperature  but  produced  particle  velocities  which 
were  inconsistent  with  the  target-temperature  Maxwell- 
Boltzmann  distribution  function,  (ii)  Within  the  next  6000 
NVT  simulation  steps,  the  frequency  of  particle-velocity  scal¬ 
ing  was  decreased  to  every  40  time  steps  while  a  Nose-Hoover 
[22]  temperature-control  algorithm  (“ thermostat ”)  was 
applied  between  the  particle -velocity  scaling  steps.  A  brief 
description  of  the  Nose-Hoover  thermostat  could  be  found  in 
our  prior  work  [23] .  (iii)  During  the  final  8000  steps,  the  tem¬ 
perature  was  controlled  using  only  the  Nose-Hoover  thermo¬ 
stat.  This  procedure  ensured  an  efficient  temperature  control 
while  yielding  an  equilibrium  state  of  the  material  (i.e.,  a 
particle-velocity  distribution  consistent  with  the  target- 
temperature  Maxwell- Boltzmann  distribution  function). 

Upon  establishing  the  thermodynamic  equilibrium  at 
5300  K,  the  target  temperature  was  reduced  by  500  K  and 
then  this  procedure  was  reapplied  at  progressively  (by  500  K) 
lower  temperatures  until  the  final  temperature  of  300  K  was 
reached.  The  total  system  equilibration  procedure  typically 
involved  simulation  times  on  the  order  of  500  ps  resulting 
in  an  average  cooling-rate  of  ~10K/ps.  A  close-up  of  the 


resulting  fused-silica  molecular-level  random  network  is 
displayed  in  Figure  2(b).  Larger  ball  sizes  are  used  in  this 
figure  to  highlight  a  pair  of  Si044-  tetrahedra  sharing  a 
common  oxygen  atom. 

2.2.  Force  Fields.  The  behavior  of  a  material  system  at  the 
molecular  level  is  governed  by  the  appropriate  force  fields 
which  describe,  in  an  approximate  manner,  the  various  inter¬ 
actions  taking  place  between  the  constituent  particles,  atoms, 
ions,  charge  groups,  and  so  forth.  In  other  words,  the  knowl¬ 
edge  of  force  fields  enables  determination  of  the  potential 
energy  of  a  system  in  a  given  configuration.  In  addition, 
gradients  of  the  force-field  functions  quantify  the  net  forces 
experienced  by  the  particles,  the  information  that  is  needed 
in  the  molecular  dynamics  simulations. 

In  general,  the  potential  energy  of  a  system  of  interacting 
particles  can  be  expressed  as  a  sum  of  the  valence  (or 
bond),  £valence,  cross-term,  Ecross_term,  and  nonbond,  £non_bond, 
interaction  energies  as 

^total  —  ^valence  ^cross-term  -^non-bond’  (-0 

The  valence  energy  generally  accounts  for  the  contribu¬ 
tion  of  valence  electrons  bonding  and  contains  the  following 
components:  (a)  a  bond  length/stretching  term;  (b)  a  two- 
bond  angle  term;  (c)  a  three-bond  dihedral/ torsion  angle 
term;  (d)  an  inversion  (or  a  four-atom  out-of-plane  interac¬ 
tion)  term;  and  (e)  the  so-called  “three-atom  Urey-Bradley 
term”  (i.e.,  the  interaction  of  two  atoms  which  are  bonded  to 
a  common  atom). 

The  cross-term  interacting  energy  accounts  for  the  cross¬ 
interactions  between  the  aforementioned  valence-energy 
components  and  includes  terms  like  (a)  stretch-stretch 
interactions  between  two  adjacent  bonds;  (b)  stretch-bend 
interactions  between  a  two -bond  angle  and  one  of  its  bonds; 
(c)  bend-bend  interactions  between  two  valence  angles  asso¬ 
ciated  with  a  common  vertex  atom;  (d)  stretch-torsion  inter¬ 
actions  between  a  dihedral  angle  and  one  of  its  end  bonds; 
(e)  stretch-torsion  interactions  between  a  dihedral  angle  and 
its  middle  bond;  (f)  bend-torsion  interactions  between  a 
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dihedral  angle  and  one  of  its  valence  angles;  and  (g)  bend- 
bend-torsion  interactions  between  a  dihedral  angle  and  its 
two  valence  angles. 

The  nonbond  interaction  term  accounts  for  the  interac¬ 
tions  between  nonbonded  atoms  and  includes  the  van  der 
Waals  energy  and  the  Coulomb  electrostatic  energy 

In  the  present  work,  the  so-called  “ COMPASS ”  (Con¬ 
densed-Phase  Optimized  Molecular  Potentials  for  Atomistic 
Simulation  Studies)  force  field  is  used  [24,  25].  This  highly 
accurate  force  field  is  of  an  ab  initio  type  since  most  of  its 
parameters  were  determined  by  matching  the  predictions 
made  by  the  ab  initio  quantum-mechanics  calculations  to 
the  condensed-matter  experimental  data.  A  summary  of  the 
COMPASS  force-field  functions  can  be  found  in  our  previous 
work  [26]. 

2.3.  Computational  Method(s).  All  the  all-atom  molecular- 
level  calculations  were  carried  out  within  the  present  work 
using  Discover  [27]  (an  atomic  simulation  program  from 
Accelrys).  Both  equilibrium  and  nonequilibrium  molecular- 
dynamics  (MD)  analyses  were  employed  in  the  present  work. 
Within  the  molecular  dynamics  method,  negative  gradient  of 
the  potential  energy  evaluated  at  the  location  of  each  atom  is 
first  used  to  compute  forces  acting  on  each  atom.  Then,  the 
associated  Newtons  equations  of  motion  (three  equations  for 
each  atom)  are  integrated  numerically  over  a  femtosecond- 
long  time  interval.  This  procedure  is  repeated  over  a  pico- 
to-nanosecond-long  simulation  time  in  order  to  determine 
the  temporal  evolution  of  the  material  molecular-level  con¬ 
figuration.  Equilibrium  MD  was  used  in  the  aforementioned 
Monte  Carlo  based  bond-switching  procedure  to  generate 
amorphous  state  from  the  crystalline  state  in  Si02.  Nonequi¬ 
librium  MD  was  used  during  the  simulation  of  the  ballistic 
impact  by  a  hard  projectile  onto  a  fused-silica  target-plate. 

Within  the  equilibrium  MD  method,  the  system  under 
consideration  is  coupled  to  an  (external)  environment  (e.g., 
constant-pressure  piston  and  constant-temperature  reser¬ 
voir)  which  ensures  that  the  system  remains  in  equilibrium 
(i.e.,  the  system  is  not  subjected  to  any  thermodynamic 
fluxes).  In  the  present  work,  NVT ,  NPT  ( P  is  pressure),  and 
NVE  ( E  is  the  total  energy)  equilibrium  MD  simulations  were 
used.  Equilibrium  MD  calculations  enable  determination  of 
the  (equilibrium)  thermodynamic  properties  of  a  material 
system  through  the  use  of  time  averages  of  the  state  variables 
sampled  along  the  calculated  system  trajectories. 

Within  nonequilibrium  MD,  the  system  is  subjected  to 
large  mechanical  and/or  thermal  perturbations  (momentum 
transfer  from  the  moving  projectile  to  the  initially  station¬ 
ary  target-plate,  in  the  present  case).  As  a  consequence, 
the  system  experiences  large  fluxes  of  its  thermodynamic 
quantities  (mass,  momentum,  and  energy,  in  the  present 
case).  Since  Discover  was  initially  designed  to  carry  out 
equilibrium  MD  simulations,  a  procedure  had  to  be  devised 
to  deactivate  “equilibration”  portions  of  this  algorithm  so 
that  nonequilibrium  MD  calculations  can  be  carried  out. 
This  procedure  was  implemented  using  a  Discover  input  file 
[25].  This  file  is  written  using  the  Basic  Tool  Command 
Language  (BTCL)  which  enabled  the  use  of  a  scripting  engine 
that  provides  very  precise  control  of  simulation  tasks  (e.g., 


specification  of  the  projectile  constant  incident  velocity  and 
deactivation  of  the  thermal- equilibration  algorithm). 

2.4.  Problem  Formulation.  The  problem  addressed  in  the 
present  work  involves  all- atom  molecular-level  computa¬ 
tional  modeling  of  the  ballistic  impact  by  a  hard  projectile 
onto  a  fused-silica  target-plate,  and  the  coordination/topo- 
logy  analysis  of  the  as-impacted  fused  silica.  The  projectile 
is  of  a  solid  right -circular  cylindrical  shape  and  impacts  the 
target -plate  normally  (i.e.,  at  a  0°  obliquity  angle  through - 
the- thickness  direction).  The  projectile  is  driven  at  a  constant 
velocity  of  1500  m/s  and  the  target-plate  is  confined  only 
along  its  bottom  rim. 

3.  Quantum-Mechanical  Analysis  of 
Glass  Devitrification 

In  order  to  rationalize  the  molecular-level  computational 
results  pertaining  to  the  ballistic  impact  of  a  projectile  onto 
the  fused-silica  target-plate,  and  the  potentially  accompa¬ 
nying  changes  in  the  glass  topology,  a  quantum-mechanical 
DFT  analysis  of  the  fused-silica  devitrification  process  is  car¬ 
ried  out  in  the  present  work.  The  main  purpose  of  this  analysis 
was  (a)  the  establishment  of  the  energy  difference  between  the 
as-received  glass  material  state  and  the  crystalline  a-quartz 
and  stishovite  Si02  states;  (b)  determination  of  the  transition 
state  energy  barriers  associated  with  the  two  (glass  — >  a- 
quartz  or  glass  — >  stishovite)  fused-silica  devitrification 
processes.  The  transition  state  is  the  material  state  along 
the  devitrification  pathway  which  is  associated  with  the 
maximum  potential  energy;  and  (c)  the  establishment  of  the 
effect  of  high  pressure  and  shear  stresses  on  the  energetics  of 
the  initial,  transition,  and  final  states. 

3.1.  Computational  Models.  As  mentioned  earlier,  coesite 
does  not  typically  form  under  ballistic-loading  conditions 
and,  hence,  will  not  be  the  subject  of  the  present  DFT 
investigation.  Consequently,  only  two  glass  devitrification 
product  states,  a-quartz  and  stishovite,  will  be  investigated. 
The  respective  computational  cells  are  shown  in  Figures  1(a) 
and  1(b).  Since  fused  silica  possesses  an  amorphous  structure, 
without  long-range  order,  a  substantially  larger  computa¬ 
tional  cell  had  to  be  used  for  this  (initial)  state  of  Si02, 
in  order  to  more  accurately  determine  its  average  potential 
energy.  An  example  of  such  a  computational  cell  is  given  in 
Figure  2(a),  the  target-plate  computational  subdomain. 

In  the  portion  of  the  analysis  dealing  with  determination 
of  the  transition  state,  the  initial  and  final  states  of  the 
material  have  to  contain  the  same  number  of  atoms  and 
species.  For  that  reason,  fused-silica  computational  cells 
containing  48  Si  and  96  O  atoms  (the  numbers  correspond  to 
a2x2x2  =  8  stack  of  the  unit  cells  depicted  in  Figure  1(a)) 
are  used  to  investigate  fused-silica  — >  a-quartz  transition 
state.  Likewise,  fused-silica  computational  cells  containing  16 
Si  and  32  O  atoms  (the  numbers  correspond  toa2x2x2  =  8 
stack  of  the  unit  cells  depicted  in  Figure  1(b))  are  used  to 
investigate  fused-silica  — >  stishovite  transition  state.  In  order 
to  account  for  the  statistical  effects  associated  with  the  extrac¬ 
tion  of  such  small  fused-silica  computational  cells  from  a 
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larger  computational  cell,  ten  small  fused- silica  computa¬ 
tional  cells  (not  shown  for  brevity)  are  generated  for  each  of 
the  two  transition  state  analyses  and  the  results  obtained  are 
averaged  out.  It  should  be  noted,  however,  that  due  to  the  fact 
that  the  initial  state  of  the  material  within  the  computational 
cells  is  amorphous,  the  cell  sizes  used  are  relatively  small, 
negatively  affecting  the  accuracy  of  the  results  obtained. 

3.2.  Computational  Method.  All  calculations  in  the  present 
work  are  carried  out  using  the  ah  initio  density-functional 
theory  code  DMol3  developed  by  Accelrys  Inc.  [28].  In 
this  code  each  electronic  wave  function  is  expanded  in  a 
localized  atom-centered  basis  set  with  each  basis  function 
defined  numerically  on  a  dense  radial  grid.  No  pseudopo¬ 
tential  approximation  is  used  for  the  near- core  electrons. 
Instead,  all-electron  calculations  are  performed  with  a  double 
numerical  polarized  (DNP)  basis  set,  the  most  complete  basis 
set  available  in  the  DMol3  code.  This  basis  set  is  equivalent 
to  the  commonly  used  analytical  6-31g**  basis  set,  a  split- 
valence  basis  set  with  polarization  functions  p  to  H  and  d  to 
C,  and  the  halogens  [29]. 

Within  the  exchange-correlation  potential  of  the  density 
functional  theory  (DFT),  the  variation  of  the  electron- charge 
density  at  the  atomic  scale  is  represented  using  the  so- 
called  Generalized  Gradient  Approximation  (GGA).  Follow¬ 
ing  Grujicic et al.  [30, 31],  the Perdew-Burke-Ernzerhof  (PBE) 
gradient-corrected  functional  [32]  is  used  in  the  present 
work.  A  standard  value  of  5.5  A  is  assigned  to  the  finite  basis- 
set  cutoff  radius. 

3.3.  Determination  of  the  Transition  States.  In  general,  the 
transition  state  lies  on  a  Minimum  Energy  Pathway  connect¬ 
ing  the  initial  and  final  states  potential  energy  minima,  where 
the  pathway  and  the  two  minima  all  reside  on  the  associated 
potential  energy  hypersurface.  The  defining  feature  of  the 
transition  state  is  that  it  is  associated  with  a  minimum 
potential  energy  in  all  the  directions  but  one  (in  which  the 
potential  energy  experiences  a  maximum).  In  the  case  of 
a  system  with  two  degrees  of  freedom,  the  transition  state 
corresponds  to  a  saddle  point,  as  depicted  in  Figure  3.  Moving 
the  system  under  consideration  from  the  transition  state 
point  (in  either  direction)  along  the  steepest-descent  path 
leads  to  the  initial/final  states. 

There  are  several  algorithms  which  are  commonly  used 
for  determination  of  the  transition  state.  The  three  most 
frequently  used  include  (a)  linear  synchronous  transit  (LST) 
[34],  used  in  the  present  work;  (b)  quadratic  synchronous 
transit  (QST)  [34];  and  (c)  nudged  elastic  band  [33]. 

The  LST  method  constructs  the  initial-state  — >  final- 
state  transition  pathway  by  connecting  each  atom  in  its 
initial  and  final  states  using  a  straight-line  pathway  and  by 
constructing  the  intermediate  configurations  by  linearly  and 
synchronously  interpolating  the  atomic  positions  along  the 
pathway  of  each  atom.  In  other  words,  the  structure  of  the 
intermediate  states  is  a  rule  of  mixtures  of  the  initial  structure 
and  the  final  structure,  with  a  weighting  factor  for  the  product 
state  /  =  (0, 1).  The  potential  energy  is  next  determined  for 
all  the  intermediate  states,  and  the  one  associated  with  the 
largest  potential  energy  is  taken  as  a  first  approximation  of  the 
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Figure  3:  The  initial-to-final  state  transition  in  a  two-degree-of- 
freedom  system  and  identification  of  the  corresponding  transition 
state  (the  saddle  point). 


system  transition  state.  The  location  of  the  transition  state  is 
further  improved  by  carrying  out  a  constrained  optimization 
of  its  first  approximation. 

4.  Results  and  Discussion 

4.1.  Validation  of  the  As-Received  Fused- Silica  Material 
State.  In  this  section,  an  attempt  is  made  to  validate  the 
fused-silica  room-temperature/ambient-pressure  structure 
obtained  through  application  of  the  previously  described 
bond-switching  and  simulated-annealing  computational 
procedures.  In  particular,  the  material  mass  density,  the 
partial  radial  distribution  functions  for  the  three  (Si-Si, 
Si-O,  and  O-O)  atomic  pairs,  and  distributions  of  the 
Si-atom  and  O-atom  coordination  numbers  are  calculated 
and  compared  with  their  experimental  counterparts. 

To  compute  the  material  mass  density,  the  partial  radial 
distribution  functions,  and  atomic  coordination  combina¬ 
tions,  NPT  equilibrium  molecular  dynamics  simulations 
were  run  at  the  room  temperature  and  the  ambient  pressure. 
Simulations  were  carried  out  for  over  10000  (0.1  fs-long)  time 
steps  while  maintaining  the  temperature  using  the  Nose- 
Hoover  thermostat  and  scaling  the  particle  velocities  every 
ten  time  steps.  Pressure  was  controlled  using  a  Berendsen 
barostat  [35].  A  brief  overview  of  this  barostat  is  provided  in 
our  prior  work  [36]. 

4.1.1.  Mass  Density.  The  average  mass  density  of  fused  silica 
computed  from  the  computational- cell  average  volume  and 
the  cell  mass  has  been  found  to  be  about  1.0%  greater  than  its 
commonly  cited  experimental  counterpart  of  2.19  g/cm3.  This 
computation/experiment  agreement  has  been  deemed  to  be 
reasonably  good. 
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Figure  4:  A  comparison  between  the  fused-silica  Si-Si,  O-O,  and  Si-O  partial  radial  distribution  functions  computed  in  (a)  the  as- 
received/initial  state,  determined  in  the  present  work,  and  (b)  the  equilibrium  initial  state  as  reported  in  Henkelman  and  Jonsson  [33]. 


4.1.2.  Radial  Distribution  Functions.  The  partial  radial  distri¬ 
bution  (often  also  referred  to  as  the  partial  pair  correlation) 
function  provides  a  measure  of  the  probability  that,  given  the 
presence  of  an  atom  of  type  a  at  the  origin  of  an  arbitrary 
reference  frame,  there  will  be  an  atom  of  type  /3  within  a 
spherical  shell  of  infinitesimal  thickness  dr  at  a  distance 
r  from  the  reference  atom.  Alternatively,  this  function  can 
be  considered  as  a  function  which  defines  a  ratio  of  the 
probability  of  finding  an  a-/3  atomic  pair  with  the  separation 
distance  r,  and  the  average  probability  of  finding  an  a-/3 
atomic  pair  at  the  same  distance.  In  amorphous  materials  like 
fused  silica,  the  partial  pair  correlation  functions  are  quite 
important  since  they  (a)  provide  an  insight  into  the  short- 
range  order  of  the  system,  (b)  can  be  used  in  the  assessment  of 
continuum-level  thermodynamic  material  properties,  and  (c) 
provide  a  way  of  validating  the  molecular-level  calculations 
since  these  quantities  can  also  be  determined  experimentally 
using  X-ray  diffraction. 

The  computed  partial  radial  distribution  functions  for 
the  as-received/ initial  state  of  fused  silica,  Figure  4(a),  are 
compared  with  their  counterparts  based  on  the  shell  model 
molecular  dynamics  calculations  [37],  Figure  4(b).  This  com¬ 
parison  reveals  that  the  present  computational  results  are 
qualitatively  similar  to  the  ones  reported  in  [38-40].  As  far 
as  the  quantitative  agreement  between  the  two  sets  of  results 
is  concerned,  it  could  be  characterized  as  being  fair  to  good. 

4.1.3.  Si-  and  O-Atom  Coordinations.  Since  fused  silica  does 
not  contain  any  network  modifiers  and  Si  is  the  only  network¬ 
forming  element,  each  Si  atom  is  expected  to  be  bonded  to 
four  O  atoms  while  each  O  atom  is  expected  to  be  bonded 
to  two  Si  atoms.  This  is  confirmed  through  postprocessing 
of  the  molecular  dynamics  results,  as  shown  in  Figures  5(a)- 
5(b)  (the  data  labeled  “initial  state”).  The  results  shown  in 


Figure  5(a)  pertain  to  the  Si-atom  coordination  while  those 
displayed  in  Figure  5(b)  refer  to  the  O-atom  coordination. 

4.2.  Analysis  of  the  Ballistic  Impact.  Temporal  evolution  of  the 
computational  domain  at  four  (0.5  ps,  1.5  ps,  2.5  ps,  and  3.5  ps) 
times  following  the  initial  contact  between  the  diamond 
solid  right- circular  cylindrical  projectile  moving  at  a  high 
velocity  and  the  fused-silica  target-plate  is  depicted  in  Figures 
6(a)-6(d),  respectively.  For  improved  clarity,  the  symbols 
for  the  target-plate  atoms  are  made  smaller.  Examination  of 
the  results  displayed  in  Figures  6(a)-6(d)  reveals  that  the 
target-plate  material  in  the  close  vicinity  of  the  projectile 
experiences  major  changes.  However,  the  nature  of  these 
changes,  that  is,  the  accompanying  alterations  in  the  material 
topology,  is  difficult  to  infer  from  the  results  displayed  in 
Figures  6(a)-6(d).  To  overcome  this  problem,  a  few  selected 
results  revealing  the  as-impacted  fused-silica  local  topology 
in  the  region  adjacent  to  the  projectile  are  presented  and 
discussed  below. 

Figure  7  shows  a  thin  slice  of  the  fused-silica  target- 
plate,  having  the  faces  parallel  to  the  x-z  plane,  centered 
on  the  hole  created  by  the  projectile  after  penetrating  the 
target  by  about  half  of  its  thickness.  Examination  of  the 
fused-silica  topology  within  this  slice  revealed  the  presence 
of  numerous  six-folded  Si  and  three-folded  O  atoms  (the 
Si  and  O  coordinations  characteristic  of  stishovite  and  not 
of  fused  silica).  In  Figure  7,  two  six-coordinated  Si  and  two 
three-coordinated  O  atoms  are  highlighted  by  assigning  a 
larger  sphere  radius  to  the  atoms  involved.  In  addition, 
six-folded  Si  atoms  and  three-folded  O  atoms  are  tagged 
with  circular  symbols.  In  addition  to  the  highlighted  atomic 
configurations,  many  more  instances  of  six-coordinated  Si 
and  three-coordinated  O  are  found  in  the  region  surrounding 
the  penetration  hole.  Since,  as  mentioned  earlier,  the  initial 
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Figure  5:  Fractional  distribution  of  the  (a)  Si-atom  and  (b)  O-atom  coordination  numbers  in  the  initial  and  the  as-impacted  states  of  fused 
silica. 
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Figure  6:  Temporal  evolution  of  the  computational  domain  at  four  times:  (a)  0.5  ps,  (b)  1.5  ps,  (c)  2.5  ps,  and  (d)  3.5  ps,  following  the  initial 
contact  between  the  diamond  solid  right- circular  cylindrical  projectile  moving  at  a  high  velocity  and  the  fused-silica  target-plate. 


state  of  fused-silica  only  contained  four-folded  Si  and  two- 
folded  O  atoms,  this  finding  suggests  that  ballistic  impact  can 
lead,  at  least  locally,  to  the  conversion  of  amorphous  Si02  into 
the  crystalline  (but  highly  deformed)  stishovite  structure. 


Additional  changes  observed  in  the  as-impacted  fused 
silica  pertain  to  the  distribution  of  the  smallest  Si-O  rings. 
To  quantify  the  size-distribution  of  the  smallest  Si-O  rings 
in  both  the  initial  and  the  as-impacted  fused-silica  states,  a 
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Figure  7:  Local  stishovite-like  topology  showing  two  six- folded  Si  and  two  three-folded  O  atoms.  For  clarity,  a  larger  sphere  radius  is  assigned 
to  the  atoms  involved.  In  addition,  six-folded  Si  atoms  and  three-folded  O  atoms  are  tagged  with  circular  symbols. 


Smallest  Si-O  ring  size  Smallest  Si-O  ring  size 

(a)  (b) 

Figure  8:  Size- distribution  function  for  the  smallest  Si-O  rings  in  the  fused-silica:  (a)  initial  state  and  (b)  as-impacted  state. 


computational  method  was  developed.  This  method  solves 
the  class  of  so-called  “shortest  path  problems”  and  is  a  simple 
modification  of  Dijkstras  algorithm  [41, 42].  The  main  modi¬ 
fications  in  this  algorithm  are  associated  with  the  fact  that,  in 
the  present  case,  the  starting  point  and  the  destination  point 
of  the  path  are  identical.  The  smallest-ring  size-distribution 
results  for  the  initial  and  as-impacted  fused-silica  states  are 
displayed  in  Figures  8 (a) -8(b),  respectively.  Examination  of 
these  results  reveals  that  ballistic  impact  alters  the  ring  struc¬ 
ture  within  fused  silica.  Specifically,  while  no  five-membered 
rings  were  present  in  the  initial  state  of  fused  silica,  such 
rings  are  found  in  the  as-impacted  state  of  the  same  material. 
Additional  changes  observed  pertain  to  the  topology  of  the 
six-membered  rings.  That  is,  while  in  the  initial  state  of  fused 
silica  these  rings  resemble  the  corresponding  rings  found  in 
cristobalite  (another  polymorph  of  Si02),  Figure  9(a),  in  the 
as-impacted  fused  silica  the  topology  of  the  six-membered 
rings  was  found  to  resemble  more  that  found  in  a-quartz, 
Figure  9(b). 

Three  partial  radial  distribution  functions  for  the  fused- 
silica  target-plate  after  the  diamond  impactor  has  penetrated 


halfway  through  the  target-plate  thickness  are  depicted  in 
Figure  10(a).  A  comparison  of  these  results  with  their  as- 
received  fused-silica  counterparts,  Figure  4(a),  reveals  that 
the  ballistic  impact  causes  distinct  changes  in  the  short-range 
order  and  atomic  coordination  within  this  material.  To  help 
rationalize  these  changes,  the  same  partial  radial  distribu¬ 
tion  functions  are  calculated  for  a-quartz,  Figure  10(b),  and 
stishovite,  Figure  10(c).  A  comparison  of  the  results  displayed 
in  Figure  4(a)  and  Figures  10 (a) -10(c)  further  confirms  that 
the  ballistic  impact-induced  changes  in  the  pair  correlation 
functions  are  associated  with  the  devitrification  of  fused 
silica,  and  the  formation  of  a-quartz  and  stishovite.  In  other 
words,  differences  in  the  partial  pair  correlation  functions 
in  the  as-impacted  fused  silica  relative  to  those  in  the  as- 
received  fused  silica  appear  to  be  caused  by  the  presence  of  Si- 
Si,  Si-O,  and  0-0  atomic  pairs  with  an  atomic  environment 
similar  to  those  found  in  stishovite  (and  a-quartz). 

To  further  demonstrate  the  conversion  of  fused  silica  to 
stishovite  as  a  result  of  the  ballistic  impact,  Si-  and  O-atom 
coordinations  in  the  fused-silica  region  surrounding  the  pen¬ 
etration  hole  are  determined.  The  results  of  this  procedure, 
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(a)  (b) 

Figure  9:  Two  views  of  the  atomic  structure  of  six-membered  Si-O  rings  in  (a)  cristobalite  and  (b)  quartz. 
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Figure  10:  Three  partial  radial  distribution  functions  for  (a)  the  fused-silica  region  adjacent  to  the  projectile,  after  the  diamond  impactor  has 
penetrated  approximately  halfway  through  the  target-plate  thickness,  (b)  a-quartz  and  (c)  stishovite. 
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Figure  11:  Relative  room-temperature  potential  energies  of  fused  silica,  a-quartz  and  stishovite,  as  a  function  of  pressure,  in  the  (a)  absence 
and  (b)  presence  of  5%  shear. 


labeled  “as-impacted  state,”  are  shown  in  Figures  5(a)-5(b). 
Examination  of  the  results  displayed  in  these  figures  reveals 
the  presence  of  six-coordinated  Si  and  three-coordinated 
O  atoms,  the  atomic  coordination  which  characterizes  the 
stishovite  crystal  structure. 

It  should  be  noted  that  all  the  results  reported  in  this 
subsection  were  obtained  under  the  conditions  of  the  pro¬ 
jectile  being  driven  at  a  constant  velocity  of  1500  m/s  during 
the  MD  simulations.  Under  such  simulation  conditions,  one 
cannot  determine  the  extent  of  momentum  transfer  from  the 
projectile  to  the  target.  In  a  separate  analysis,  the  projectile 
was  assigned  an  initial  velocity  of  1500  m/s  rather  than  being 
driven  at  this  velocity.  In  this  simulation,  it  was  found  that 
about  60%  of  the  incident  linear  momentum  of  the  projectile 
is  transferred  to  the  target.  In  other  words,  the  residual 
velocity  of  the  projectile  after  passing  through  the  target  is 
approximately  600  m/s. 

4.3.  Relative  Stability  of  Fused  Silica ,  a-Quartz  and  Stishovite. 
To  further  help  rationalize  the  changes  in  the  partial  pair 
correlation  functions  induced  by  the  ballistic  impact,  the 
quantum-mechanical  DFT  calculation  results  pertaining  to 
the  relative  stability  of  fused  silica,  a-quartz  and  stishovite, 
as  well  as  of  the  energy  barriers  associated  with  the  fused- 
silica  — >  a-quartz  and  fused-silica  — >  stishovite  transition 
states,  are  presented  and  discussed  in  this  section.  The  relative 
room-temperature  potential  energies  of  the  fused  silica,  a- 
quartz  and  stishovite  as  a  function  of  pressure  are  presented 
in  Figure  11(a).  Examination  of  the  results  displayed  in  this 
figure  reveals  that,  at  the  ambient  pressure,  a-quartz  is  the 


most  stable  form  of  Si02,  followed  by  fused  silica  and  then 
by  stishovite.  At  a  pressure  of  50  GPa,  on  the  other  hand,  the 
relative  positions  of  a-quartz  and  fused  silica,  with  respect  to 
thermodynamic  stability,  have  been  exchanged. 

In  Figure  10(b),  the  relative  room-temperature  potential 
energies  of  the  fused  silica,  a-quartz  and  stishovite  as  a 
function  of  pressure,  in  the  presence  of  5%  shear,  are 
presented.  A  comparison  of  the  results  displayed  in  Figures 
11(a) -11(b)  reveals  that  the  application  of  5%  shear  strain 
did  not  change  the  relative  stability  of  fused  silica,  a-quartz 
and  stishovite  at  the  ambient  pressure.  In  sharp  contrast,  at 
50  GPa  pressure,  the  application  of  5%  shear  has  been  found 
to  make  the  stishovite  the  most  stable  phase,  followed  by  a- 
quartz  and  then  fused  silica.  This  finding  suggests  that,  in 
the  presence  of  shear,  fused  silica  is  more  likely  to  undergo 
a  devitrification  conversion  into  stishovite  and,  to  a  lower 
extent,  into  a-quartz.  This  finding  is  important  since  the 
fused-silica  region  surrounding  the  projectile  is  generally 
subjected  to  high  shear,  in  addition  to  high  pressures.  This 
explains  why,  in  Figure  7,  this  region  was  found  to  undergo 
an  extensive  fused-silica  — >  stishovite  conversion. 

The  likelihood  for  the  aforementioned  devitrification 
processes  is  affected  not  only  by  the  relative  stabilities  of  fused 
silica,  a-quartz  and  stishovite,  but  also  by  the  size  of  the 
energy  barrier  associated  with  the  respective  transition  state. 
Variations  in  the  fused-silica  — »  a-quartz  and  fused-silica  — > 
stishovite  transition  state  energy  barriers  with  pressure,  in  the 
absence  and  the  presence  of  5%  shear,  are  shown,  respectively, 
in  Figures  12(a) -12(b).  Examination  of  the  results  displayed 
in  these  figures,  at  a  pressure  of  50  GPa,  reveals  that  the 
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Figure  12:  Variations  in  the  fused-silica  — >  a-quartz  and  fused-silica  — >■  stishovite  transition  state  energy  barriers  with  pressure,  in  the  (a) 
absence  and  (b)  presence  of  5%  shear. 


energy  barrier  for  fused-silica  devitrification,  in  the  absence 
of  shear,  is  slightly  higher  for  the  fused-silica  — >  stishovite 
conversion.  On  the  other  hand,  the  energy  barrier  for  fused- 
silica  devitrification,  in  the  presence  of  shear,  is  slightly  higher 
for  the  fused-silica  — >  a-quartz  conversion.  These  findings 
suggest  that,  in  the  presence  of  shear,  and  at  pressures  as  high 
as  50  GPa,  fused  silica  is  more  likely  to  convert  to  stishovite 
than  a-quartz. 

Figure  13  shows  the  conversion  of  an  initially  amorphous 
Si02  structure  into  a  stishovite -like  structure  under  the 
influence  of  high  pressure  and  shear.  The  transition  state 
associated  with  this  conversion  is  also  shown  in  this  figure. 
For  clarity,  only  a  small  Si02  region  consisting  of  three  Si  and 
fifteen  O  atoms  is  displayed  in  Figure  13.  Due  to  such  a  small 
size  of  the  region,  some  O  atoms  appear  to  be  nonbonded. 
These  O  atoms  are  bonded  to  Si  atoms,  but  the  Si  atoms  that 
they  are  bonded  to  are  not  shown  in  this  figure.  The  Si  atoms 
in  question  are  also  bonded  to  some  of  the  bonded  O  atoms 
displayed  in  Figure  13.  Examination  of  this  figure  reveals 
that,  as  expected,  the  fused-silica  state  contains  only  fourfold 
coordinated  Si  and  twofold  coordinated  O  atoms,  while  the 
stishovite  state  contains  sixfold  coordinated  Si  and  threefold 
coordinated  O  atoms.  The  Si-  and  O-atom  coordination  in 
the  transition  state  involves  fourfold  and  fivefold  coordinated 
Si  and  twofold  and  threefold  coordinated  O  atoms.  To  help 
with  the  interpretation  of  the  results  displayed  in  Figure  13, 
the  same  three  Si  atoms  appearing  in  the  three  configurations 
are  denoted  using  labels  “Si  1,”  “Si  2,”  and  “Si  3.”  It  should  be 
recalled  that  the  results  presented  in  Figures  12(a) -12(b)  reveal 
that  the  maximum  energy  associated  with  the  transition  state 
is  lowered  and,  thus,  the  conversion  of  the  fused-silica  — > 
stishovite  becomes  feasible,  under  high  pressure  and  shear. 


Figure  13:  Conversion  of  an  initially  amorphous  Si02  structure  into 
a  stishovite-like  structure  under  the  influence  of  high  pressure  and 
shear.  The  associated  transition  state  is  also  shown. 


5.  Summary  and  Conclusions 

Based  on  the  results  obtained  in  the  present  work,  the 
following  summary  remarks  and  main  conclusions  can  be 
drawn. 

(1)  By  determining  the  Si-Si,  Si-O,  and  0-0  partial  radial 
distribution  functions,  and  the  coordination  numbers  for  Si 
and  O  atoms,  before  and  after  a  ballistic  impact  experienced 
by  a  fused-silica  target-plate,  potential  devitrification  of  the 
initial  material  and  its  conversion  into  a-quartz  and  stishovite 
are  investigated. 
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(2)  The  investigation  is  carried  out  computationally 
through  the  use  of  all-atom  nonequilibrium  molecular- 
dynamics  simulations. 

(3)  To  rationalize  the  results  obtained,  quantum-mech¬ 
anical  density  functional  theory  computational  analyses  are 
employed.  These  analyses  helped  to  establish  the  effect  of  high 
pressure  and  shear  on  the  relative  stability  of  the  amorphous 
fused  silica,  crystalline  a-quartz  and  crystalline  stishovite, 
and  on  the  energy  barriers  associated  with  the  fused-silica  — > 
a-quartz  and  fused-silica  — >  stishovite  conversion  reactions. 

(4)  The  results  obtained  demonstrated  that  under  pres¬ 
sures  on  the  order  of  50  GPa,  and  in  the  presence  of  shear, 
stishovite  becomes  the  most  stable  form  of  Si02  and  the 
energy  barrier  associated  with  the  fused-silica  — >  stishovite 
conversion  is  the  smallest.  Since  the  fused-silica  target- 
plate  regions  beneath  and  surrounding  the  hard  projectile 
experience  high  pressures  and  shear,  formation  of  stishovite 
and,  to  a  lower  extent,  a-quartz  from  fused  silica,  as  observed 
in  the  all- atom  molecular-level  computational  analyses  of  the 
ballistic-impact  problem,  is  justified. 
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